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THE SUBGRID PROBLEM OF THE THERMAL CONVECTION 
IN THE EARTH'S LIQUID CORE 

M. Reshetnyak^^^ g StefFen^ 
Abstract 

The problem of the turbulent thermal convection in the Earth's liquid core is considered. Following 
assumptions on decreasing of the spatial scales due to the rapid rotation, we propose the subgrid 
model of the eddy diffusivity, which is used in the large-scale model. This approach makes it possible 
to model realistic regimes with small Ekman and Rossby numbers {E ^ 10^^'', Ro ^ 10^^) and 
a sufficiently large Rayleigh number Ra ~ 10^^. The obtained estimate of the averaged kinetic 
energy is comparable with observations. The model includes rotation of the solid core due to the 
viscous torque. 
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1 Introduction 

Convection in the hquid core of the Earth, caused by the radiactive heating and compositional processes 

is the subject of numerous researches, usually concerned with the geomagnetic field generation, also. 

The last few decades saw a fascinating development in this area 0. Based on the MHD large-scale 

equations, numerical simulations can reproduce different geomagnetic and geophysical phenomena: 

various properties of the geomagnetic field (e.g., its reversals and spectrum), eastward rotation of the 

inner core of the Earth as well as the realistic ratio of the kinetic and magnetic energies |3] • 

However, the wide range of spatial and temporal scales make the direct numerical simulations 

(DNS) very cumbersome. The difficulty is caused by the small values of the transport coefficients: for 

instance the kinematic viscosity of the liquid core is = 10~^m^s~^, and the thermal diffusivity: 

= 10^^ m^s~^ (here the subscript ^ corresponds to the molecular values). This gives estimates 
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of the molecular Reynolds and Peclet numbers of i?e*^ = ^^^^5^ ~ 10^ and Pe*^ = ^^^^^ ~ 10^, where 
V^d = 0.2° year~^ is the west drift velocity and L = 3 • 10^ m is the scale of the liquid core 
which corresponds to the regime of the highly developed turbulence. In the case of the Kolmogorov's 
turbulence 3D DNS, simulations require ~ = 10^^ grid nodes Attempts to use the exact 

values of these parameters on the coarse grid lead to the numerical instabilities. The first intuitive 
models in geodynamo theory which suppressed instabilities at the small scales, e.g., the model of 
hyperdiffusivity [S] , gave rise to new questions concerned with interpretation of the results obtained P| . 
The more consequential way is an application of the semiempirical models of turbulence . Usually, 
these models are based on assumptions on the cascade transfer similar to Kolmogorov's, which give 
descriptions of the average effect of the small-scale field fluctuations onto the large-scale flow in terms 
of the eddy diffusivity. The recent studies of the subgrid and complex models jHlEl of the thermal 
convection and dynamo problems in the rotating sphere revealed the principal possibility of describing 
the small-scale fluctuations in the turbulence with the desired Reynolds and Peclet numbers much like 
Kolmogorov's model. 

These models work up to the regime of moderate rotation speed. Further increase of the Coriolis 
force can reduce the total kinetic energy and even suppress convection at all. From linear analysis it 
follows that the critical Rayleigh number depends on the Ekman number like ~ E~^^^ Even 
though the molecular estimate of the Rayleigh number gives huge numbers R^^ ~ 10^^ [HIIS], this 
value is only 5 • 10^ times larger then the critical value R^ 0. Due to the rapid rotation of the Earth, 
the situation in the liquid core is more complicated and assumptions on similarity of the spectral 
characteristics of the fields must be checked very carefully. We show that the direct applications of 
the traditional models of the turbulence, based on the mix-length assumptions, lead to results that 
differ from the observations by orders of magnitude. The cause of such disagreement lies in the daily 
rapid rotation of the Earth, which gives rise to new characteristic spatial scales in the core As a 
result, the energy distribution in the spectrum changes, which makes the application of the Prandtl- 
Kolmogorov's approach to the eddy diffusivity estimate difficult. Convection at these new scales plays 
the crucial role in the energy balance of the whole system and changes the estimate of the total energy 
by orders of magnitude. Even a simple account of these effects leads to essential change of the rate of 
the energy dissipation and thus to a better agreement of LSS models with the observations. 

In the section ^ we introduce the large-scale equations of the thermal convection and consider 
the Prandtl-Kolmogorov's assumptions on the eddy diffusion. In the section |21 we recall the basics of 
convection in a rapidly rotating body and estimate the subgrid diffusion. Afterwards, this estimate is 
used in the large-scale model, section^ The discussion of results is in sectional 
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2 The large-scale equations 



The problem of the thermal convection in the Earth's core can be reduced to the problem in the 
spherical shell. Let the surface of the sphere, radius ro (in the spherical system of coordinates (r, 6, 93)), 
rotating with angular velocity around the z-axis. This sphere contains a concentric solid inner sphere, 
radius r^, and the outer spherical layer (rj < r < tq) is filled with an incompressible liquid (V = 0). 
The inner sphere is allowed to rotate freely around the z-axis due to viscous torque. Convection in 
the Boussinesq approximation in the outer sphere is described by the Navier-Stokes equation and by 
the heat flux equation. Choosing L = tq as the unit of length, velocity V, time t and pressure p can 
be measured in units of k^^ /L, L'^/k^'^ and 2Q,pK^^ , respectively. Then, the governing equations can 
be written in the form 

'dV 



^^+(V.V)V 



-Vp + -UxY + R^Trlr + E^^V- S, (1) 



^ + (V.V)(r + To)=V-(VT), (2) 



where Iz is the unit vector in z-direction, S is the rate of the strain tensor and T the temperature 

"i/r-i 
1-n 



fluctuations relatively to the imposed proflle Tq = ^]^^^ . The molecular Rossby i?*^ , Ekman E^^ and 



Rayleigh numbers appear in the equations 

pM _ k'^' rpM _ v'"\ r>M _ agpSTL (o\ 

where a is the coefficient of thermal expansion, g the gravity acceleration and 6T is a temperature 
unit {5T ^ 10~^K, see [2]). It should be mentioned that the Rayleigh number for non-rotating bodies 
is usually given in the form = agST / u'^ and that R^^ = E^^R^^. 

The solid inner sphere is allowed to rotate freely around the z-axis due to viscous torque. The 
dimensionless momentum equation for the angular velocity u) of the inner sphere (0 < r < r^) has the 
form 



Rol^ = riE^' i> Vlr=r, Sin ^ dS, (4) 



'dt 

s 

where / is the moment of inertia of the inner sphere S and 5^^ is a component of the strain tensor in the 
spherical system of coordinates . Equations (^^21 B]) are accompanied by the non-penetrating and 
no-slip boundary condition for velocity V and zero temperature fluctuations at the shell boundaries. 

The system (fTl l2l4[) was successfully studied in the regimes of the laminar convection using different 
numerical approaches jl6l . However, these regimes are still very far from the desired estimates for 
the Earth's hquid core R^ = W'^, E^^ = 10" and R^^ = 10^^ 2 . Attempts to approach to these 
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parameters using DNS caused numerical instabilities and required application of turbulencs models 
[S]. However, even the direct usage of the known models of turbulence is not trivial. 

To support this point we offer a simple estimate of the eddy diffusivity, based on the most popular 
mix-length model of the turbulence. Following the Prandtl-Kolmogorov hypotheses, the eddy diffusion 
at the scale / can be estimated as z^^ = (e/^)^/^, where e = v'^ /I is the rate of energy dissipation and v 
is a velocity at the scale I. Even the largest estimate, based on the main scale / = L and the west drift 
velocity ^ = 3 • 10~'^ms~^, gives = 2 - lO^m^s"^, giving an Ekman number of order = 2 • 10~^. 
The more realistic estimate with Vi = 6V ~ ^(7;)^^^ and the usual grid scale of Z ~ 3 • 10"^L gives 
iJ^ = 15 and E'^ = 10~^. On the other hand, this estimate of E'^ would require resolution of about 

~ 2-10^ columns which need use of the most powerful modern computers. All this means that 
this estimate of iJ^ will not provide the smooth field behaviour of fields assumed in the Kolmogorov's 
turbulence, when E'^ > 1. Thus, the traditional methods underestimate the eddy diffusion u^. 

Such situation corresponds to the case, where the classical ideas on the direct cascade of energy 
from the main scale L to the dissipative scale are violated and additional information on e at the 
dissipative scale is needed. As we see below in the section |31 it appears that in the case of the rapid 
rotating body the energy in the spectrum is shifted to the small scales, those which DNS cannot 
resolve even at the onset of convection. This is the reason, why any attempts to estimate v'^ in the 
turbulent regime at scales compared with the grid resolution, lead to the non-selfconsistent behaviour 
of the turbulent model. 

The way out of such difficulties is make proper assumptions on the spectral properties of the 
solution in the range of the high wave numbers. 

3 The model of columns 

The origin of the problem can be seen from the analysis of the linerized system (|1I2|) at the onset of 
convection in the limit of small Rossby and Ekman numbers. As it was shown in 11 already (see 
also recent paper ^\), at the onset of convection the structure of the flow tends to develop columns 
along z-direction, such that d/d'^ ~ 0{E~^/^), d/ds ~ 0{E-^'^), d/dz ~ 0(1), when E = Ro 
Linerization of the system (^EJ leads to the balance of the Archemedean and viscous terms in the 
Navier-Stokes equation: i?a T ~ E~^/^V . The balance of the convective and viscous terms in the 
heat-flux equation gives V ^ E^'^^^T, from which follows the estimate of the critical Rayleigh number 
R^^ ~ £^-1/3. (For convenience we omitted index .) Such, at the onset of convection for system HH2() . 
the flow is anisotropical with the smallest scale Ie ~ E^^^L, defined by the balance of the Coriolis and 
viscous forces. Note that the scale Ie ~ 10~^ is beyond the level of DNS. If this asymptotic is correct, 
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the critical Rayleigh number in the Earth's core is R"^' ~ 10^ P]- As we show below, the predicted 
column-like form of the flow is very important for estimates of the subgrid dissipation in the liquid 
core. 

The main assumption is, that even in the turbulent regime believed to be in the Earth's liquid 
core, the flow tends to elongated structures with the smallest scale Ie, predicted by the linear analysis. 
It is from this scale ideas of the direct cascade of energy are applicable. To simplify the problem, we 
estimate the isotropical eddy diffusion, based on the scale {Ie)- In particular, instead of the estimate 
of velocity gradient at the subgrid scale I: V' '-^ SVjl^ we use y E~^/^6V, where 6V ~ 0.31^ 
is the average variation of velocity at the scale /. In this case the estimate of the eddy diffusion 
gives ~ PV ~ 5 • Iff^m^ s~^ and E'^ ~ 5 • 10~^. This estimate of the turbulent Ekman number 
corresponds to ~ 10 columns which can be resolved in the large-scale models with the desired 
accuracy. To demonstrate these arguments, we propose a simulations of the system ^ with the 
given eddy diffusion i/^ estimated as above. 



4 Turbulent model. Results of calculations. 

Equations (^12 EJ) are solved using the control- volume method (Simple algorithm) on the staggered 
grid (rir, uq, n^) = (45,45,64). This method is based on the finite-difference approximation and 
demonstrates very high numerical stability for the regimes with strong convection^. For ease of 
calculation, we renormalize equations 0J) using turbulent diffusion units, so that instead of 

the K = Im^ s""*^ was used. Then, the dimensionless parameters are: = ~ ^ " I0~^' 

= = 10"'^. We consider three regimes with turbulent Rayleigh numbers i?J = " ^"n^^ ~ ^^^^ 
10^, 10^ (see the time evolution of the kinetic energy Ek in Fig. The corresponding Reynolds 
numbers averaged over the shell volume Re = -^y^IEk are 3 • 10^, 6 • 10^ and 2 • 10^°, c. f. with 
the molecular Reynolds number for the Earth's core based on the west drift velocity R^'^^^^ ~ 10^. 

Characteristic snapshots of the large-scale velocity r, 6, (/^-components are presented in Fig. |21 The 
observed curls in r, (/^-projections corresponds to the columns parallel to z-axis. These columns may 
drift in the (/9-direction. In its turn, the non-zero viscous gradient u-^ (^) causes rotation of the 
inner core, (see evolution of the angular velocity of the inner core lo in Fig^. Here the positive value 
of uj corresponds to the eastward direction, known to occur in the Earth [201 • emphasize that 
these maps are a product of averaging of the small-scale {Ie ~ 0{E^'^/^) = 10 ^) structures. So far, 
the micro-scale Reynolds number re at the scale Ie is still larger then unity, and the inertial spectrum 
for the scales smaller then Ie exists. An estimate of = with / = E^^^^^ and v = O.IV gives 
^See also some special questions of the control- volume method for the full dynamo problem in the sphere in |19|. 
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Tg ~ 10^. This spectrum has two parts with the transition point defined by the balance of the inertial 
and Coriohs terms: Iq ~ RqV. The turbulence in the range of Z^; < /j^ is influenced by rotation and 
the kinetic energy spectrum is Ei ~ P |21j . For the scales smaller then Iq up to the dissipative scale 
Id = Re "^^^ the Kolmogorov's spectrum Ei ~ 1^1'^ reappears. 

Summarizing the obtained results we conclude, that based on the realistic values of the Rossby 
and Rayleigh numbers and on assumptions on the spectrum of the flow in the liquid core we obtained 
a value of the kinetic energy Ek comparable with the observations. Having in mind that the velocity 
field and the eddy diffusion are connected in our model, we consider this agreement to be worth to 
note. 

5 Conclusions 

We propose the scenario of the turbulent thermal convection in the rapid rotating body, when the 
Coriolis force shifts the system to the origin of the small scales already at the onset of convection, and 
show that further increasing of the intensity of the heat sources leads to a turbulent regime, which is 
still far from the Kolmogorov's case. It appears that predictions of the linear analysis at the onset of 
convection are applicable to the eddy diffusion estimate in the regime of the fully developed turbulence. 
Though the original problem is highly anisotropical, the "isotropical" estimate of the eddy diffusion 
gives a kinetic energy of the system comparable with the observations. Note that introduction of 
the magnetic field will not change the problem in principal, because at the scales Ie considered the 
corresponding micro-scale magnetic Reynolds number is already *C 1 and the magnetic field decays 
due to the Ohmic dissipation process. On the other hand, it is not yet clear how the west drift velocity 
relates to the flow at the scales Ie and different interpretations of observations can exist. This question 
requires the solution of the full dynamo problem. 
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Figure 1: Evolution of the angular velocity of the liquid core u and averaged over the volume kinetic 
energy Ek for i?o = 4 • 10"^, E = 10"=^; Rl = IQ^ - thick line, (2) - Rl = lO'^ - thin line, (3) - 
R^ = IQS - dashed line. 
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Figure 2: The snapshots of the velocity field components {vr, v^, v^) (from left to right) for the equa- 
torial sections (the left half of the plane): (-700, 1200), (-4700, 4100), (-1200, 1800) and merid- 
ional sections for axi-symmetrical parts of the fields (the right half): (—2400, 1800), (—3800, 3100), 
(—400, 1000). Numbers in round brackets indicate ranges. 
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